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ABSTRACT 

The availability and proper utilization of fluid properties is of fundamental importance in the 
process of mathematical modeling of propulsion systems. Real fluid properties provide the bridge 
between the realm of pure analytics and empirical reality. The two most common approaches used to 
formulate thermodynamic properties of pure substances are fundamental (or characteristic) equations of 
state (Helmholtz and Gibbs functions) and a piecemeal approach that is described, for example, in 
Adebiyi and Russell (1992). This paper neither presents a different method to formulate thermodynamic 
properties of pure substances nor validates the aforementioned approaches. Rather its purpose is to 
present a method to be used to facilitate the accurate interpretation of fluid thermodynamic property data 
generated by existing property packages. There are two parts to this paper. The first part of the paper 
shows how efficient and usable property tables were generated, with the minimum number of data points, 
using an aerospace industry standard property package (based on fundamental equations of state 
approach). The second part describes an innovative interpolation technique that has been developed to 
properly obtain thermodynamic properties near the saturated liquid and saturated vapor lines. 

INTRODUCTION 

For many years, computer-based mathematical models have been used to analyze 
thermodynamic systems, and it is no doubt an effective way to save time and money for doing parametric 
studies and/or what-if scenarios. The quality of a model depends on many factors such as assumptions, 
governing equations, boundary condition definition, and of course thermodynamic properties. Having 
good fluid property data not only improves the model’s accuracy, but it also makes the model easier to 
debug. For example, when the modeling results are in question relative to, say, hardware testing results, 
it is necessary to determine if there are any errors in the model. If the fluid property routines are not 
trusted, then they may be contributing factors to the errors; however, if fluid property routines have been 
verified, then they are no longer suspects, and the modeler can be confident that the errors are 
elsewhere. 

The fundamental (or characteristic) equations of state and the “piecemeal” method are the two 
most common approaches used to formulate thermodynamic properties of pure substances. The 
fundamental equations of state are derived from a combination of the first and second laws of 
thermodynamics along with the Helmholtz and the Gibbs functions. Using these laws and functions, the 
complete list of thermodynamic properties can then be obtained by differentiating the fundamental 
equations of state. On the other hand, the thermodynamic properties equations for the “piecemeal” 
approach, which are developed by using the specific heat and pressure, specific volume, and 
temperature data for each phase of a substance, are typically the integrals of the pertinent functions. 

Ideally, a simple set of equations would be integrated into a code (just like conservation of mass, 
momentum, and energy) and used to obtain every thermodynamic property for all phases of a given 
substance, but, unfortunately, that is not the case. So the next best thing is to use existing property 
packages, which contain complex equations with many coefficients, to get the fluid property data 
necessary for use within the mathematical models. However, it is very time consuming to use standard 
property packages within a computer simulation; therefore, property tables are used to save 
computational time. With that in mind, there are two parts to this paper. The first part of the paper shows 
how efficient and usable property tables were generated, with the minimum number of data points, using 
an aerospace industry standard property package (based on fundamental equations of state approach). 
The second part describes an innovative interpolation technique that has been developed to properly 
obtain thermodynamic properties near the saturated liquid and saturated vapor lines. 



RESULTS AND DISCUSSION 


PROPERTY PACKAGES 


There are a number of property packages commercially available. NIST, MIPROPS, and 
GASPAK are a few examples of those packages and were used in this paper for comparison testing. As 
a demonstration of the techniques discussed in this paper, selected sets of data points for oxygen were 
chosen to represent every part of an entire fluid region. Non-saturated properties were generated using 
pressure (P) and temperature (T) as inputs, while saturated properties were obtained using pressure and 
quality (X). Table 1 compares density (p) and constant pressure specific heat (C p ) from the three property 
packages listed above, and as expected, the thermodynamic properties for those data points are very 
similar. Thus, for this example, it is a matter of personal preference as to which property package to use 
to generate property tables. The authors have decided to use GASPAK due to its relative ease of 
implementation. 



NIST 

MIPROPS 

GASPAK 

P = 220 psia; X = 0 
Saturated Liquid 

p (Ibm/ft 3 ) 

57.6702 

57.6710 

57.6687 

C D (Btu/lbm-°R) 

0.4951 

0.4950 

0.4951 

P = 220 psia; X = 1 
Saturated Vapor 

p (Ibm/ft^) 

3.6691 

3.6760 

3.6688 

C D (Btu/lbm-°R) 

0.3548 

0.3500 

0.3548 

P = 520 psia; T = 200 °R 
Subcooled Liquid 

p (Ibm/ft 3 ) 

64.9515 

65.0110 

64.9514 

Co (Btu/lbm-°R) 

0.4238 

0.4220 

0.4238 

P = 70 psia; T = 500 °R 
Superheated Vapor 

p (Ibm/ft^) 

0.4192 

0.4200 

0.4192 

C D (Btu/lbm-°R) 

0.2208 

0.2210 

0.2208 

P = 5000 psia; T = 700 °R 
Supercritical 

p (Ibm/ft 3 ) 

19.6892 

19.7160 

19.6889 

C D (Btu/lbm-°R) 

0.2755 

0.2760 

0.2755 


Table 1 . Property Data for Different Property Packages 


PROPERTY TABLES 


The first step in creating computationally efficient fluid property tables involves the evaluation of 
the property state space. Real fluid property tables were generated using the GASPAK property package 
(based on fundamental equations of state approach) by inputting pressures and enthalpies. The 
fundamental goal of this effort was to minimize the number of pressure-enthalpy values necessary to 
generate property tables that fully described the fluid. To accomplish this, a computer program was 
created which indirectly allows the user to determine how many pressure-enthalpy pairs are necessary in 
order to achieve a prescribed interpolation error. Thus, the desired accuracy of the results across the 
region of fluid properties considered is used to determine the expanse of the tables generated. 

First, the program evaluates fluid properties along the saturated liquid and vapor lines to 
determine the necessary number of pressure points and their values so that the interpolation error along 
these lines would be within the specified error. This is illustrated in Figure 1 (a). Then the fluid property 
data is divided into two regions, sub-critical and supercritical, for further evaluation. After obtaining the 
pressure values for saturated lines, the program then determines pressure-enthalpy values for the 
supercritical region. This is a more involved process than what have been done for the saturated lines, 
but more or less, the methodology remains functionally similar. The program starts out with the lowest 
enthalpy value and determines the necessary number of pressure points for the supercritical region so 
that the interpolation error along that enthalpy line would be within the specified error. Similarly, the 
program starts out with the critical pressure (lowest pressure point for the supercritical region) and 
determines the necessary number of enthalpy points required to fulfill the specified interpolation error 
limits. The process is repeated up and across the entire supercritical region until the grid is generated as 
shown in Figure 1 (b). This same technique could be used for the sub-critical region, and one would find 
that the sub-critical enthalpies are almost identical to the supercritical enthalpies. So for the sub-critical 


region, it is acceptable to use the same pressure values as the saturated lines, and likewise, the same 
enthalpy values as the supercritical region. The entire grid is shown in Figure 1 (c). The aforementioned 
program could be used to determine the grid for each fluid property at a time or to determine the grid for 
multiple fluid properties at the same time. In other words, one could use the program to obtain the grid for 
density and then use the program again to obtain the grid for specific heat, and most likely, the two grids 
would not be the same. On the other hand, one could use the program to determine only one grid that 
would satisfy the specified interpolation error limits for both density and specific heat. Figure 2 
demonstrates what the grids would look like for oxygen temperature and density, and as expected, the 
grids get finer in the vicinity of the critical point. 

INTERPOLATION METHODS 

After having evaluated the state space for a given fluid property, the data is compiled into tables. 
Then, a means of utilizing these tables in an accurate manner across the full state space must be applied. 
Two tables were generated for each thermodynamics property. The first table (Property Table 1) only 
includes saturated vapor and saturated liquid properties, while the second table (Property Table 2) 
contains everything else but saturated properties. For the supercritical region and much of the space of 
the sub-critical region, the table-reading method used is straightforward double interpolation of Table 2. 
For points within the dome (or saturation envelope), data along the saturation lines (Property Table 1) are 
used along with fluid quality to determine composite fluid properties. However, for a small part of the sub- 
critical region, a special interpolation technique is required for cases where one or more data points within 
Property Table 2 falls inside the saturation dome. Figure 3 shows fourteen possible cases near the 
saturated lines where the special interpolation method is required in order to retrieve reasonable property 
data from the generated tables. 

To see how the special interpolation technique works, let us consider saturated liquid Case 2 as 
defined in Figure 3 and as shown in Figure 4 within the state space. The top two data points within the 
interpolation rectangle (enclosed by four unshaded circular dots - from Property Table 2) lie inside the 
saturation envelope; therefore, they will not be used in the interpolation process. Instead, only the bottom 
two data points are used along with two other data points from the saturated liquid line (unshaded square 
dots - from Property Table 1). By using the input pressure (P) and data from Property Table 1 , the 
saturated properties are then calculated and shown as a shaded square dot on the saturated liquid line. 

In the same way, the other properties (shaded circular dot) are computed by using the same input 
pressure and data from Property Table 2. Finally, the properties of interest are determined by 
interpolating between the newfound properties using the input enthalpy (H). This methodology is also 
used for the other thirteen cases defined in Figure 3; however, some of the other cases are not as 
straightforward as the example case. 

Now let us look at another case, saturated vapor Case 6 as defined in Figure 3 and as shown in 
Figure 5 within the state space, where it is not as straightforward as the saturated liquid Case 2. Again, 
the interpolation rectangle is formed by four unshaded circular dots (from Property Table 2), and the 
bottom left data point will not be used in the interpolation process because it lies inside the saturated 
envelope. First, it is necessary to determine the pressure (Pcross) at which the saturated vapor line 
crosses the lower enthalpy line (H to ) of the interpolation rectangle and then to see if the input pressure (P) 
is higher than or lower than Pcross- So within saturated vapor Case 6, there are two sub-cases. One case 
is where P is less than Pcross, and the other case is where P is greater than or equal to Pcross- This 
example case falls into the latter sub-case. Next by using Pcross and data from Property Table 1 , the 
saturated properties are calculated and shown as a shaded square dot on the saturated vapor line. Now 
the property for H| 0 (shaded triangular dot) is computed by interpolating between the newfound saturated 
properties and the bottom right data point from Property Table 2 using the input pressure. Similarly using 
the same input pressure, the property for H w (shaded circular dot) is obtained by interpolating between 
the top two data points from Property Table 2, and finally the properties of interest are determined by 
interpolating between the H| 0 and H M properties using the input enthalpy (H). 




Log P 
(c) 


Figure 1 . Pressure-Enthalpy Grid 
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Figure 2. Property Data for Oxygen Temperature and Density 
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Figure 3. Possible Special Interpolation Cases Near Saturated Lines 
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Figure 4. Saturated Liquid Case 2 
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Figure 5. Saturated Vapor Case 6 
PROPERTY TABLES VERSUS PROPERTY PACKAGES 


As stated earlier, the main reason for using property tables instead of property packages within a 
computer simulation is to save computational time, but how much time is actually saved by using property 
tables? To answer this question, a simple model was built where there is only one control volume (node), 
and the thermodynamic properties for that volume were determined 50,000 times by using the same 
pressure and enthalpy each time. When property tables were used, it only took about 13 seconds to get 
the oxygen properties; however, it took almost 25 minutes (24 minutes and 54 seconds) when GASPAK 
was used to perform the same task. One might be surprised by how long it took when using the property 
package, but remember that GASPAK was called within the computer simulation. There is no doubt that 
it would take a lot less time if GASPAK is used outside of the computer simulation, but the basis for 
building property tables is not to use them outside of a computer simulation, but to use them within a 
computer simulation. 

Now let us consider a mathematical model for a hypothetical facility used for component testing 
and see how long it would take for the model to run. Let us assume that the model is comprised of 25 
control volumes. Fifteen of those volumes have only one fluid, either a fuel or an oxidizer, and the other 
ten volumes have two fluids, a fuel and an inert fluid or an oxidizer and an inert fluid. First, let us consider 
the one-fluid volumes. In order to get property data for a one-fluid volume, a computer code that would 
be used in this example needs to call a property package (or property tables) three times, so the total 
number of property call for fifteen one-fluid volumes is 45 times. For a two-fluid volume, the computer 
code needs to call a property package six times, three times per fluid. Since there are two fluids, it is 
assumed that the fluids would mix together and become a homogeneous mixture, and in order to get 
reasonable properties for the mixture, it is assumed that on average it would take about three iterations. 
So for a two-fluid volume, it would require a total of 18 property calls, and consequently, the total number 
of property call for ten two-fluid volumes is 180 times. Hence, the total number of property call for this 
twenty-five-volume example is 225 times. Now let us assume that it would take, on average, three 
iterations for the overall system model to converge for each time step, so the total number of property call 




per converged time step is 675 times. Finally, if the simulation is run for five seconds with a time step of 

0.001 second then the total number of property call is 3,375,000 times. As noted earlier for every 50,000 
property calls, it would take about 1 3 seconds when using property tables and about 25 minutes when 
using property package. So for this example, it would take about 1 5 minutes versus 28 hours to finish a 
five seconds simulation, and that is a lot of time saving. This example is illustrated in Table 2, and it 
represents a fairly typical mathematical modeling situation. 


1 

One-Fluid Volume 

Two-Fluid Volume 


Number of property call for one volume 

3 

6 

One 

Homogeneous mixture - Number of iteration 

1 

3 

Time 

Total number of property call for one volume 

3 

18 

Step 

Number of volume in the model 

15 

10 



45 

180 


Total number of property call for the model 

225 


Number of iteration for each time step to converge 

3 


Total number of property call for each time step 

675 


Number of time step 

(5 seconds simulation, 0.001 second time step) 

5000 


Total number of property call for the simulation 

3,375,000 



Property Tables 



Time required for every 50,000 property calls 

13 sec. 

24 min. and 54 sec. 


Total simulation time 

-15 min. 

-28 hrs. 


Table 2. Twenty-Five-Volume Model 

SUMMARY AND CONCLUSIONS 

There are two parts to this paper. The first part of the paper shows how efficient and usable 
property tables were generated, with the minimum number of data points, using the GASPAK property 
package. The numbers of data points in the tables were indirectly determined from the specified 
acceptable error between interpolated values from the tables and values obtained from the property 
package. The second part of the paper describes how thermodynamics properties were retrieved from 
the tables. For the most part, the interpolation method used is straightforward use of property tables. 
However, for the fourteen cases near the saturation lines, the special interpolation method is required in 
order to obtain reasonable property data from the generated tables. As a result, a new interpolation 
routine was developed to handle both sub-critical and supercritical regions. 
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